Mark D. Scheuerell
Fish Ecology Division, Northwest Fisheries Science Center, National Marine Fisheries Service, National Oceanic and Atmospheric Administration, Seattle, WA USA, mark.scheuerell@noaa.gov

0.0.1 Version

This is version 0.17.08.01.


1 Initialization

library(readr)
library(reshape2)
library(MARSS)

2 Data

For these analyses we will use time series of counts of adult spring/summer Chinook salmon from the Snake River basin in Oregon and Idaho. This group of fish forms an Evolutionarily Significant Unit, which was listed as threatened under the US Endangered Species Act in 1992.

The data were compiled largely through the efforts of Tom Cooney (Northwest Fisheries Science Center, National Marine Fisheries Service, National Oceanic and Atmospheric Administration, Portland, OR USA, tom.cooney@noaa.gov).

2.1 Assets

The assets data file has the following columns of information:

  • pop = abbreviated name of the population
  • code = 5-letter abbreviation for MPG (2 letters) & population (3 letters)
  • mpg = name of 1 of 5 Major Population Groups
  • ha = estimated area (hectares) of available spawning habitat
  • category = hatchery supplementation category (“impact” or “control”)
  • brood_yr = year when spawning occurred
  • nS = total number of spawning fish
  • phos = proportion of hatchery-origin spawners
rawdata <- read.csv("https://raw.githubusercontent.com/mdscheuerell/CAPM/master/data/srss_chin.csv")
head(rawdata)
##      pop  code                   mpg      ha category brood_yr   nS phos
## 1 CathCr GRCAT Grande Ronde / Imnaha 29.1522   impact     1955  122    0
## 2 CathCr GRCAT Grande Ronde / Imnaha 29.1522   impact     1956 2606    0
## 3 CathCr GRCAT Grande Ronde / Imnaha 29.1522   impact     1957 1754    0
## 4 CathCr GRCAT Grande Ronde / Imnaha 29.1522   impact     1958  486    0
## 5 CathCr GRCAT Grande Ronde / Imnaha 29.1522   impact     1959  712    0
## 6 CathCr GRCAT Grande Ronde / Imnaha 29.1522   impact     1960 3161    0

The first thing to do is remove hatchery-origin spawners from the counts, and then convert the counts to log-density based on the estimated spawning area for each population. In a couple of years, incomplete censuses with small sample sizes led to estimated fractions of wild fish equal to zero. I will treat those as NA.

dat <- rawdata[,c("mpg", "code", "pop", "category", "brood_yr")]
## calculate log-density
dat$Sw <- round(log(rawdata$nS*(1-rawdata$phos)/rawdata$ha), 2)
## replace possible -Inf with NA
dat$Sw[is.infinite(dat$Sw)] <- NA

For easier bookkeeping, I’ll also switch the Imnaha’s MPG over to the Grande Ronde.

## change Imnaha name for easier sorting later
dat[,"code"] <- sub("IRMAI", "GR_IR", dat[,"code"], fixed=TRUE)
## MPG abbreviations
mpg_ID <- c("GR", "MF", "SF", "SR")

Only a few of the populations have observations going all the way back to the 1950s, so I’ll trim the data accordingly.

## use generation time of 4 years
gt <- 4
## first brood year to consider
yr_first <- 1960
## last brood year to consider
yr_last <- 2011
## years of interest
t_index <- seq(yr_first, yr_last-gt)
## length of time period
TT <- length(t_index)
## subset data
dat <- subset(dat, brood_yr>=yr_first & brood_yr<=yr_last)

I will use the MARSS package to fit the CAPM model, which requires the data to be formatted as an [N x T] matrix.

## reshape assets
assets <- dcast(dat, brood_yr ~ code, value.var="Sw")[,-1]
## compute differences
assets <- t(apply(assets, 2, diff, lag=gt))
matplot(t(assets), type="b")

## number of popns
NN <- dim(assets)[1]
## names of popns
names_pop <- rownames(assets)

Finally, here are some nicer, longer names for the various populations.

# set up naming & ordering scheme
longnames <- c(
  "Wenaha",
  "Grande Ronde",
  "Catherine",
  "Minam",
  "Lostine",
  "Imnaha",
  "SF Salmon",
  "Secesh",
  "EFSF Salmon",
  "Big",
  "Sulphur",
  "Bear Valley",
  "Marsh",
  "Loon",
  "Camas",
  "Yankee Fork",
  "Valley",
  "UM Salmon",
  "LM Salmon",
  "EF Salmon",
  "Lemhi")
name_ord <- c(6,3,5,4,2,1,12,10,15,14,13,11,9,7,8,20,21,19,18,17,16)
names_mat <- data.frame(ord=name_ord, ID=names_pop)
names_mat <- names_mat[order(names_mat$ord),]
names_mat$long <- longnames

2.2 Market index

We used the monthly PDO data provided by the University of Washington’s Joint Institute for the Study of the Atmosphere and Ocean (JISAO), which are available here. We begin by downloading the raw PDO data and viewing the metadata.

## raw PDO data from UW
PDO_raw <- readLines("http://research.jisao.washington.edu/pdo/PDO.latest")
## line with data headers
hdr_pdo <- which(lapply(PDO_raw,grep,pattern="YEAR")==1, arr.ind=TRUE)
## print PDO metadata
print(PDO_raw[seq(hdr_pdo-1)],quote=FALSE)
##  [1]                                                                                                                                   
##  [2] PDO INDEX                                                                                                                         
##  [3]                                                                                                                                   
##  [4] If the columns of the table appear without formatting on your browser, use http://research.jisao.washington.edu/pdo/PDO.latest.txt
##  [5]                                                                                                                                   
##  [6] Updated standardized values for the PDO index, derived as the                                                                     
##  [7] leading PC of monthly SST anomalies in the North Pacific Ocean,                                                                   
##  [8] poleward of 20N. The monthly mean global average SST anomalies                                                                    
##  [9] are removed to separate this pattern of variability from any                                                                      
## [10] "global warming" signal that may be present in the data.                                                                          
## [11]                                                                                                                                   
## [12]                                                                                                                                   
## [13] For more details, see:                                                                                                            
## [14]                                                                                                                                   
## [15]  Zhang, Y., J.M. Wallace, D.S. Battisti, 1997:                                                                                    
## [16]      ENSO-like interdecadal variability: 1900-93. J. Climate, 10, 1004-1020.                                                      
## [17]                                                                                                                                   
## [18]  Mantua, N.J. and S.R. Hare, Y. Zhang, J.M. Wallace, and R.C. Francis,1997:                                                       
## [19]      A Pacific interdecadal climate oscillation with impacts on salmon                                                            
## [20]      production. Bulletin of the American Meteorological Society, 78,                                                             
## [21]      pp. 1069-1079.                                                                                                               
## [22]      (available via the internet at url:                                                                                          
## [23]          http://www.atmos.washington.edu/~mantua/abst.PDO.html)                                                                   
## [24]                                                                                                                                   
## [25] Data sources for this index are:                                                                                                  
## [26]  UKMO Historical SST data set for 1900-81;                                                                                        
## [27]  Reynold's Optimally Interpolated SST (V1) for January 1982-Dec 2001)                                                             
## [28] *** OI SST Version 2 (V2) beginning January 2002 -                                                                                
## [29]                                                                                                                                   
## [30] missing value flag is -9999                                                                                                       
## [31]                                                                                                                                   
## [32]

Next, we will extract the actual PDO indices for the years of interest and inspect the file contents.

## PDO data for years of interest
dat_PDO <- read.table("http://research.jisao.washington.edu/pdo/PDO.latest",
                      header=FALSE, stringsAsFactors=FALSE,
                      skip=hdr_pdo + yr_first - 1900 + 1 + 2, nrows=TT)
colnames(dat_PDO) <- unlist(strsplit(tolower(PDO_raw[hdr_pdo]), split="\\s+"))
dat_PDO
##      year   jan   feb   mar   apr   may   jun   jul   aug   sep   oct   nov   dec
## 1    1962 -1.29 -1.15 -1.42 -0.80 -1.22 -1.62 -1.46 -0.48 -1.58 -1.55 -0.37 -0.96
## 2    1963 -0.33 -0.16 -0.54 -0.41 -0.65 -0.88 -1.00 -1.03  0.45 -0.52 -2.08 -1.08
## 3    1964  0.01 -0.21 -0.87 -1.03 -1.91 -0.32 -0.51 -1.03 -0.68 -0.37 -0.80 -1.52
## 4    1965 -1.24 -1.16  0.04  0.62 -0.66 -0.80 -0.47  0.20  0.59 -0.36 -0.59  0.06
## 5    1966 -0.82 -0.03 -1.29  0.06 -0.53  0.16  0.26 -0.35 -0.33 -1.17 -1.15 -0.32
## 6    1967 -0.20 -0.18 -1.20 -0.89 -1.24 -1.16 -0.89 -1.24 -0.72 -0.64 -0.05 -0.40
## 7    1968 -0.95 -0.40 -0.31 -1.03 -0.53 -0.35  0.53  0.19  0.06 -0.34 -0.44 -1.27
## 8    1969 -1.26 -0.95 -0.50 -0.44 -0.20  0.89  0.10 -0.81 -0.66  1.12  0.15  1.38
## 9    1970  0.61  0.43  1.33  0.43 -0.49  0.06 -0.68 -1.63 -1.67 -1.39 -0.80 -0.97
## 10   1971 -1.90 -1.74 -1.68 -1.59 -1.55 -1.55 -2.20 -0.15  0.21 -0.22 -1.25 -1.87
## 11   1972 -1.99 -1.83 -2.09 -1.65 -1.57 -1.87 -0.83  0.25  0.17  0.11  0.57 -0.33
## 12   1973 -0.46 -0.61 -0.50 -0.69 -0.76 -0.97 -0.57 -1.14 -0.51 -0.87 -1.81 -0.76
## 13   1974 -1.22 -1.65 -0.90 -0.52 -0.28 -0.31 -0.08  0.27  0.44 -0.10  0.43 -0.12
## 14   1975 -0.84 -0.71 -0.51 -1.30 -1.02 -1.16 -0.40 -1.07 -1.23 -1.29 -2.08 -1.61
## 15   1976 -1.14 -1.85 -0.96 -0.89 -0.68 -0.67  0.61  1.28  0.82  1.11  1.25  1.22
## 16   1977  1.65  1.11  0.72  0.30  0.31  0.42  0.19  0.64 -0.55 -0.61 -0.72 -0.69
## 17   1978  0.34  1.45  1.34  1.29  0.90  0.15 -1.24 -0.56 -0.44  0.10 -0.07 -0.43
## 18   1979 -0.58 -1.33  0.30  0.89  1.09  0.17  0.84  0.52  1.00  1.06  0.48 -0.42
## 19   1980 -0.11  1.32  1.09  1.49  1.20 -0.22  0.23  0.51  0.10  1.35  0.37 -0.10
## 20   1981  0.59  1.46  0.99  1.45  1.75  1.69  0.84  0.18  0.42  0.18  0.80  0.67
## 21   1982  0.34  0.20  0.19 -0.19 -0.58 -0.78  0.58  0.39  0.84  0.37 -0.25  0.26
## 22   1983  0.56  1.14  2.11  1.87  1.80  2.36  3.51  1.85  0.91  0.96  1.02  1.69
## 23   1984  1.50  1.21  1.77  1.52  1.30  0.18 -0.18 -0.03  0.67  0.58  0.71  0.82
## 24   1985  1.27  0.94  0.57  0.19  0.00  0.18  1.07  0.81  0.44  0.29 -0.75  0.38
## 25   1986  1.12  1.61  2.18  1.55  1.16  0.89  1.38  0.22  0.22  1.00  1.77  1.77
## 26   1987  1.88  1.75  2.10  2.16  1.85  0.73  2.01  2.83  2.44  1.36  1.47  1.27
## 27   1988  0.93  1.24  1.42  0.94  1.20  0.74  0.64  0.19 -0.37 -0.10 -0.02 -0.43
## 28   1989 -0.95 -1.02 -0.83 -0.32  0.47  0.36  0.83  0.09  0.05 -0.12 -0.50 -0.21
## 29   1990 -0.30 -0.65 -0.62  0.27  0.44  0.44  0.27  0.11  0.38 -0.69 -1.69 -2.23
## 30   1991 -2.02 -1.19 -0.74 -1.01 -0.51 -1.47 -0.10  0.36  0.65  0.49  0.42  0.09
## 31   1992  0.05  0.31  0.67  0.75  1.54  1.26  1.90  1.44  0.83  0.93  0.93  0.53
## 32   1993  0.05  0.19  0.76  1.21  2.13  2.34  2.35  2.69  1.56  1.41  1.24  1.07
## 33   1994  1.21  0.59  0.80  1.05  1.23  0.46  0.06 -0.79 -1.36 -1.32 -1.96 -1.79
## 34   1995 -0.49  0.46  0.75  0.83  1.46  1.27  1.71  0.21  1.16  0.47 -0.28  0.16
## 35   1996  0.59  0.75  1.01  1.46  2.18  1.10  0.77 -0.14  0.24 -0.33  0.09 -0.03
## 36   1997  0.23  0.28  0.65  1.05  1.83  2.76  2.35  2.79  2.19  1.61  1.12  0.67
## 37   1998  0.83  1.56  2.01  1.27  0.70  0.40 -0.04 -0.22 -1.21 -1.39 -0.52 -0.44
## 38   1999 -0.32 -0.66 -0.33 -0.41 -0.68 -1.30 -0.66 -0.96 -1.53 -2.23 -2.05 -1.63
## 39   2000 -2.00 -0.83  0.29  0.35 -0.05 -0.44 -0.66 -1.19 -1.24 -1.30 -0.53  0.52
## 40   2001  0.60  0.29  0.45 -0.31 -0.30 -0.47 -1.31 -0.77 -1.37 -1.37 -1.26 -0.93
## 41 2002**  0.27 -0.64 -0.43 -0.32 -0.63 -0.35 -0.31  0.60  0.43  0.42  1.51  2.10
## 42 2003**  2.09  1.75  1.51  1.18  0.89  0.68  0.96  0.88  0.01  0.83  0.52  0.33
## 43 2004**  0.43  0.48  0.61  0.57  0.88  0.04  0.44  0.85  0.75 -0.11 -0.63 -0.17
## 44 2005**  0.44  0.81  1.36  1.03  1.86  1.17  0.66  0.25 -0.46 -1.32 -1.50  0.20
## 45 2006**  1.03  0.66  0.05  0.40  0.48  1.04  0.35 -0.65 -0.94 -0.05 -0.22  0.14
## 46 2007**  0.01  0.04 -0.36  0.16 -0.10  0.09  0.78  0.50 -0.36 -1.45 -1.08 -0.58
## 47 2008** -1.00 -0.77 -0.71 -1.52 -1.37 -1.34 -1.67 -1.70 -1.55 -1.76 -1.25 -0.87
## 48 2009** -1.40 -1.55 -1.59 -1.65 -0.88 -0.31 -0.53  0.09  0.52  0.27 -0.40  0.08

Notice that values in the year column in the PDO data file contains an ** superscript from 2002 onward to indicate the switch to version 2 of the sea surface temperature (SST) data. We will fix those as well.

## fix 'year' col
dat_PDO$year <- seq(1:TT) + yr_first - 1
## mean PDO from Apr-Oct indexed by brood year
PDO_ts <- apply(dat_PDO[,-1],1,mean)
dat_PDO <- data.frame(year=dat_PDO$year, PDO=PDO_ts)
## market index
market <- round(dat_PDO$PDO, 3)
plot.ts(market)

2.3 Risk-free index

For these purposes, the risk free index is zero for all years. The idea is that any population with a standardized return of zero is replacing itself and therefore not declining.

free <- rep(0, TT)

3 CAPM in MARSS

3.1 v1: a static; b unconstrained

## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL

## number of regr params = level + slope
mm <- 2

## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
  
  ## print loop ID
  print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
  ## get popns
  assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
  ## number of assets
  nn <- dim(assets_mpg)[1]
  
  ## Process eqn
  ## B is Identity
  BB <- "identity"
  ## U is col vec of 0's
  UU <- "zero"
  ## Q will be block diagonal
  QQ <- matrix(list(0),mm*nn,mm*nn)
  ## upper L block for alpha
  aa <- 1:nn
  ## lower R block for beta
  bb <- aa + nn
  ## equal var-cov
  # QQ[bb,bb] <- "beta.cov"
  # diag(QQ[bb,bb]) <- "beta.var"
  ## diagonal & unequal
  diag(QQ[bb,bb]) <- paste0("beta",seq(nn))
  for(j in 1:(nn-1)) {
    for(k in (j+1):nn) {
      QQ[j+nn,k+nn] <- paste0("beta",j,k)
      QQ[k+nn,j+nn] <- paste0("beta",j,k)
    }
  }
  
  ## Observation eqn
  ## Z is array of regr variables (ie, 1's and market returns)
  ZZ <- array(NA, c(nn,mm*nn,TT))
  for(t in 1:TT) {
    ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
  }
  ## A is col vec of 0's
  AA <- "zero"
  ## assume same obs var & cov
  ##RR <- matrix(list(0),nn,nn)
  ##RR[aa,aa] <- "obs.cov"
  ##diag(RR[aa,aa]) <- "obs.var"
  RR <- "diagonal and equal"
  
  ## Initial states
  ##x0 <- rbind(matrix("pi.a",nn,1), matrix("pi.b",nn,1))
  ##V0 <- 100*diag(mm*nn)
  ##x0 <- rbind(matrix(0,nn,1), matrix(1,nn,1))
  ##V0 <- matrix(list(0),mm*nn,mm*nn)
  ##diag(V0[aa,aa]) <- "var.a"
  ##diag(V0[bb,bb]) <- "var.b"
  
  ## starting values for regr parameters
  ini_list <- list(x0=matrix(1, mm*nn, 1))
  
  ## list of model matrices & vectors
  mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
  ##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
  
  ## list of control values
  con_list <- list(safe=TRUE, allow.degen=FALSE,
                   conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
  
  ## fit multivariate DLM
  mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
  
  ## get list of Kalman filter output
  kf_out = MARSSkfss(mod_res[[i]])
  
  ## ts of forecasted betas
  betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
  
  ## ts of SE of betas; nxT matrix
  betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
  
  ## ts of alphas
  alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
  ## ts of SE{alphas}
  alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
  
  ## ts of betas
  betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
  ## ts of SE{betas}
  betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
  
}
## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1301 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1301 iterations. 
## Log-likelihood: -372.3136 
## AIC: 812.6272   AICc: 822.3414   
##  
##          Estimate
## R.diag     0.5000
## Q.beta1    1.9085
## Q.beta12   1.7668
## Q.beta13   2.5422
## Q.beta14   2.1624
## Q.beta15   2.2498
## Q.beta16   2.3169
## Q.beta2    1.8029
## Q.beta23   2.2418
## Q.beta24   1.8395
## Q.beta25   1.9753
## Q.beta26   2.0533
## Q.beta3    3.4610
## Q.beta34   2.9889
## Q.beta35   3.0685
## Q.beta36   3.1475
## Q.beta4    2.6078
## Q.beta45   2.6534
## Q.beta46   2.7143
## Q.beta5    2.7211
## Q.beta56   2.7902
## Q.beta6    2.8631
## x0.X1      0.0506
## x0.X2      0.1166
## x0.X3      0.2123
## x0.X4      0.1170
## x0.X5      0.0219
## x0.X6      0.0690
## x0.X7      0.2993
## x0.X8      1.5459
## x0.X9     -0.0543
## x0.X10    -0.2095
## x0.X11    -0.3033
## x0.X12     0.1886
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 2224 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2224 iterations. 
## Log-likelihood: -397.1267 
## AIC: 862.2534   AICc: 872.4243   
##  
##          Estimate
## R.diag     0.4909
## Q.beta1   12.9265
## Q.beta12  14.1924
## Q.beta13  11.1113
## Q.beta14  11.7245
## Q.beta15  12.3104
## Q.beta16  11.0612
## Q.beta2   15.6232
## Q.beta23  12.3696
## Q.beta24  13.0052
## Q.beta25  13.5118
## Q.beta26  12.0437
## Q.beta3   10.3303
## Q.beta34  10.6704
## Q.beta35  10.5630
## Q.beta36   9.0572
## Q.beta4   11.0871
## Q.beta45  11.1514
## Q.beta46   9.6881
## Q.beta5   11.7242
## Q.beta56  10.5450
## Q.beta6    9.7275
## x0.X1      0.0847
## x0.X2      0.0837
## x0.X3     -0.1976
## x0.X4     -0.0232
## x0.X5      0.0689
## x0.X6      0.0163
## x0.X7     -0.0795
## x0.X8      0.2604
## x0.X9     -0.3900
## x0.X10     0.0398
## x0.X11    -0.1894
## x0.X12    -0.1460
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 655 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 655 iterations. 
## Log-likelihood: -145.5705 
## AIC: 317.1409   AICc: 320.1245   
##  
##          Estimate
## R.diag    0.26835
## Q.beta1   1.24348
## Q.beta12  0.98983
## Q.beta13  0.91831
## Q.beta2   0.78821
## Q.beta23  0.72955
## Q.beta3   0.68857
## x0.X1    -0.00702
## x0.X2    -0.07964
## x0.X3     0.13335
## x0.X4     0.67648
## x0.X5     0.75019
## x0.X6     0.93047
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1560 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1560 iterations. 
## Log-likelihood: -372.5814 
## AIC: 813.1629   AICc: 822.9571   
##  
##          Estimate
## R.diag     0.4367
## Q.beta1    8.7974
## Q.beta12   7.5370
## Q.beta13   6.6868
## Q.beta14   8.8338
## Q.beta15   9.2384
## Q.beta16  11.6973
## Q.beta2    6.4641
## Q.beta23   5.7336
## Q.beta24   7.5611
## Q.beta25   7.8883
## Q.beta26  10.0038
## Q.beta3    5.0860
## Q.beta34   6.7095
## Q.beta35   7.0034
## Q.beta36   8.8786
## Q.beta4    8.8781
## Q.beta45   9.3050
## Q.beta46  11.7648
## Q.beta5    9.8061
## Q.beta56  12.3538
## Q.beta6   15.6003
## x0.X1      0.2527
## x0.X2      0.1038
## x0.X3      0.0971
## x0.X4      0.2392
## x0.X5      0.1699
## x0.X6      0.2439
## x0.X7      0.2062
## x0.X8      0.2658
## x0.X9      0.4366
## x0.X10    -0.1898
## x0.X11    -0.3480
## x0.X12    -0.4036
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop

3.1.1 Plots of results

betas <- betas[,names_mat$ID]

# total num of ts
NN <- dim(betas)[2]

# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
             rainbow(3, start=0.2, end=0.45, v=0.7),
             rainbow(6, start=0.8, end=0.95, v=0.85),
             rainbow(6, start=0.55, end=0.7))

# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7

# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))

# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)

# draw plots
for(i in 1:NN) {
    plot(t_index, betas[,i], type="n",
         ylim=yext, xaxt="n", yaxt="n",
         xlab="", ylab="", main=names_mat[i,"long"])
    rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
    abline(h=0, lty="dashed")
    box()
    lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
    lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
    lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
    #text(xin, yin, names.pop[i], adj=c(0,1))   
    if(i<=nr | nc==1) {
#       axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
#       axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
        axis(side=2, las=1, tick=TRUE)
        }
    if(i%%nr==0) {
        axis(side=1, at=t_index[(t_index %% 10)==0])
        }
    }
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)

3.2 v1: a static; b equalvarcov

## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL

## number of regr params = level + slope
mm <- 2

## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
  
  ## print loop ID
  print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
  ## get popns
  assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
  ## number of assets
  nn <- dim(assets_mpg)[1]
  
  ## Process eqn
  ## B is Identity
  BB <- "identity"
  ## U is col vec of 0's
  UU <- "zero"
  ## Q will be block diagonal
  QQ <- matrix(list(0),mm*nn,mm*nn)
  ## upper L block for alpha
  aa <- 1:nn
  ## lower R block for beta
  bb <- aa + nn
  ## equal var-cov
  QQ[bb,bb] <- "beta.cov"
  diag(QQ[bb,bb]) <- "beta.var"

  ## Observation eqn
  ## Z is array of regr variables (ie, 1's and market returns)
  ZZ <- array(NA, c(nn,mm*nn,TT))
  for(t in 1:TT) {
    ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
  }
  ## A is col vec of 0's
  AA <- "zero"
  ## assume same obs var & cov
  ##RR <- matrix(list(0),nn,nn)
  ##RR[aa,aa] <- "obs.cov"
  ##diag(RR[aa,aa]) <- "obs.var"
  RR <- "diagonal and equal"
  
  ## Initial states
  ##x0 <- rbind(matrix("pi.a",nn,1), matrix("pi.b",nn,1))
  ##V0 <- 100*diag(mm*nn)
  ##x0 <- rbind(matrix(0,nn,1), matrix(1,nn,1))
  ##V0 <- matrix(list(0),mm*nn,mm*nn)
  ##diag(V0[aa,aa]) <- "var.a"
  ##diag(V0[bb,bb]) <- "var.b"
  
  ## starting values for regr parameters
  ini_list <- list(x0=matrix(1, mm*nn, 1))
  
  ## list of model matrices & vectors
  mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
  ##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
  
  ## list of control values
  con_list <- list(safe=TRUE, allow.degen=FALSE,
                   conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
  
  ## fit multivariate DLM
  mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
  
  ## get list of Kalman filter output
  kf_out = MARSSkfss(mod_res[[i]])
  
  ## ts of forecasted betas
  betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
  
  ## ts of SE of betas; nxT matrix
  betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
  
  ## ts of alphas
  alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
  ## ts of SE{alphas}
  alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
  
  ## ts of betas
  betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
  ## ts of SE{betas}
  betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
  
}
## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1380 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1380 iterations. 
## Log-likelihood: -378.8426 
## AIC: 787.6852   AICc: 789.5034   
##  
##            Estimate
## R.diag       0.6063
## Q.beta.var   2.0867
## Q.beta.cov   2.0866
## x0.X1        0.0302
## x0.X2        0.0503
## x0.X3        0.1612
## x0.X4        0.1015
## x0.X5       -0.0105
## x0.X6        0.0184
## x0.X7        0.1434
## x0.X8        0.3911
## x0.X9        0.2313
## x0.X10       0.5211
## x0.X11       0.0687
## x0.X12       0.4557
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 1390 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1390 iterations. 
## Log-likelihood: -404.0779 
## AIC: 838.1558   AICc: 840.053   
##  
##            Estimate
## R.diag       0.6716
## Q.beta.var   9.1373
## Q.beta.cov   9.1372
## x0.X1        0.0181
## x0.X2        0.0015
## x0.X3       -0.1194
## x0.X4       -0.0641
## x0.X5        0.0122
## x0.X6       -0.0411
## x0.X7       -0.0973
## x0.X8        0.0320
## x0.X9       -0.2713
## x0.X10      -0.1796
## x0.X11      -0.2088
## x0.X12      -0.0799
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 1029 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1029 iterations. 
## Log-likelihood: -147.174 
## AIC: 312.348   AICc: 313.7765   
##  
##             Estimate
## R.diag      0.303214
## Q.beta.var  0.753774
## Q.beta.cov  0.753725
## x0.X1       0.000199
## x0.X2      -0.035161
## x0.X3       0.082217
## x0.X4       0.634399
## x0.X5       0.898669
## x0.X6       0.885139
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1394 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1394 iterations. 
## Log-likelihood: -387.9918 
## AIC: 805.9837   AICc: 807.8157   
##  
##            Estimate
## R.diag       0.5493
## Q.beta.var   7.4219
## Q.beta.cov   7.4218
## x0.X1        0.2113
## x0.X2        0.0932
## x0.X3        0.1140
## x0.X4        0.2188
## x0.X5        0.1759
## x0.X6        0.0990
## x0.X7        0.0423
## x0.X8       -0.0261
## x0.X9        0.2312
## x0.X10      -0.1561
## x0.X11       0.1660
## x0.X12      -0.2863
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop

3.2.1 Plots of results

betas <- betas[,names_mat$ID]

# total num of ts
NN <- dim(betas)[2]

# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
             rainbow(3, start=0.2, end=0.45, v=0.7),
             rainbow(6, start=0.8, end=0.95, v=0.85),
             rainbow(6, start=0.55, end=0.7))

# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7

# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))

# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)

# draw plots
for(i in 1:NN) {
    plot(t_index, betas[,i], type="n",
         ylim=yext, xaxt="n", yaxt="n",
         xlab="", ylab="", main=names_mat[i,"long"])
    rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
    abline(h=0, lty="dashed")
    box()
    lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
    lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
    lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
    #text(xin, yin, names.pop[i], adj=c(0,1))   
    if(i<=nr | nc==1) {
#       axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
#       axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
        axis(side=2, las=1, tick=TRUE)
        }
    if(i%%nr==0) {
        axis(side=1, at=t_index[(t_index %% 10)==0])
        }
    }
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)

3.3 v1: a static; b diag & unequal

## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL

## number of regr params = level + slope
mm <- 2

## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
  
  ## print loop ID
  print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
  ## get popns
  assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
  ## number of assets
  nn <- dim(assets_mpg)[1]
  
  ## Process eqn
  ## B is Identity
  BB <- "identity"
  ## U is col vec of 0's
  UU <- "zero"
  ## Q will be block diagonal
  QQ <- matrix(list(0),mm*nn,mm*nn)
  ## upper L block for alpha
  aa <- 1:nn
  ## lower R block for beta
  bb <- aa + nn
  ## diag & unequal
  diag(QQ[bb,bb]) <- paste0("beta",seq(nn))

  ## Observation eqn
  ## Z is array of regr variables (ie, 1's and market returns)
  ZZ <- array(NA, c(nn,mm*nn,TT))
  for(t in 1:TT) {
    ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
  }
  ## A is col vec of 0's
  AA <- "zero"
  ## assume same obs var & cov
  ##RR <- matrix(list(0),nn,nn)
  ##RR[aa,aa] <- "obs.cov"
  ##diag(RR[aa,aa]) <- "obs.var"
  RR <- "diagonal and equal"
  
  ## Initial states
  ##x0 <- rbind(matrix("pi.a",nn,1), matrix("pi.b",nn,1))
  ##V0 <- 100*diag(mm*nn)
  ##x0 <- rbind(matrix(0,nn,1), matrix(1,nn,1))
  ##V0 <- matrix(list(0),mm*nn,mm*nn)
  ##diag(V0[aa,aa]) <- "var.a"
  ##diag(V0[bb,bb]) <- "var.b"
  
  ## starting values for regr parameters
  ini_list <- list(x0=matrix(1, mm*nn, 1))
  
  ## list of model matrices & vectors
  mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
  ##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
  
  ## list of control values
  con_list <- list(safe=TRUE, allow.degen=FALSE,
                   conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
  
  ## fit multivariate DLM
  mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
  
  ## get list of Kalman filter output
  kf_out = MARSSkfss(mod_res[[i]])
  
  ## ts of forecasted betas
  betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
  
  ## ts of SE of betas; nxT matrix
  betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
  
  ## ts of alphas
  alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
  ## ts of SE{alphas}
  alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
  
  ## ts of betas
  betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
  ## ts of SE{betas}
  betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
  
}
## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1276 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1276 iterations. 
## Log-likelihood: -420.0602 
## AIC: 878.1204   AICc: 881.0434   
##  
##          Estimate
## R.diag   1.014417
## Q.beta1  0.016075
## Q.beta2  0.300899
## Q.beta3  0.096780
## Q.beta4  0.000323
## Q.beta5  0.028867
## Q.beta6  0.012729
## x0.X1   -0.058006
## x0.X2    0.038213
## x0.X3   -0.033504
## x0.X4   -0.018423
## x0.X5   -0.070038
## x0.X6   -0.054581
## x0.X7   -0.080212
## x0.X8    1.108690
## x0.X9   -0.737083
## x0.X10  -0.014390
## x0.X11  -0.107951
## x0.X12   0.222093
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 2156 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2156 iterations. 
## Log-likelihood: -459.2719 
## AIC: 956.5437   AICc: 959.5959   
##  
##          Estimate
## R.diag   1.676423
## Q.beta1  0.000370
## Q.beta2  0.054952
## Q.beta3  0.000115
## Q.beta4  0.012216
## Q.beta5  0.017199
## Q.beta6  0.050396
## x0.X1   -0.086090
## x0.X2   -0.039976
## x0.X3   -0.233259
## x0.X4   -0.059126
## x0.X5   -0.150412
## x0.X6   -0.322697
## x0.X7   -0.192926
## x0.X8    0.684116
## x0.X9   -0.455433
## x0.X10   0.079412
## x0.X11  -0.152422
## x0.X12  -0.205091
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 108 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 108 iterations. 
## Log-likelihood: -172.9275 
## AIC: 365.8549   AICc: 367.6149   
##  
##         Estimate
## R.diag    0.4606
## Q.beta1   0.4191
## Q.beta2   0.2020
## Q.beta3   0.1118
## x0.X1    -0.0528
## x0.X2    -0.0472
## x0.X3     0.1730
## x0.X4     0.3972
## x0.X5     1.0640
## x0.X6     0.7844
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 2756 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2756 iterations. 
## Log-likelihood: -479.2527 
## AIC: 996.5053   AICc: 999.4511   
##  
##          Estimate
## R.diag   1.64e+00
## Q.beta1  2.11e-04
## Q.beta2  9.32e-05
## Q.beta3  1.01e-04
## Q.beta4  5.99e-05
## Q.beta5  1.37e-01
## Q.beta6  9.92e-01
## x0.X1   -1.09e-01
## x0.X2   -2.27e-01
## x0.X3   -1.74e-01
## x0.X4   -1.01e-01
## x0.X5   -2.90e-01
## x0.X6   -3.54e-01
## x0.X7   -2.20e-01
## x0.X8   -2.92e-01
## x0.X9   -1.36e-02
## x0.X10  -4.21e-01
## x0.X11  -3.34e-01
## x0.X12  -1.21e+00
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop

3.3.1 Plots of results

betas <- betas[,names_mat$ID]

# total num of ts
NN <- dim(betas)[2]

# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
             rainbow(3, start=0.2, end=0.45, v=0.7),
             rainbow(6, start=0.8, end=0.95, v=0.85),
             rainbow(6, start=0.55, end=0.7))

# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7

# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))

# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)

# draw plots
for(i in 1:NN) {
    plot(t_index, betas[,i], type="n",
         ylim=yext, xaxt="n", yaxt="n",
         xlab="", ylab="", main=names_mat[i,"long"])
    rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
    abline(h=0, lty="dashed")
    box()
    lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
    lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
    lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
    #text(xin, yin, names.pop[i], adj=c(0,1))   
    if(i<=nr | nc==1) {
#       axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
#       axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
        axis(side=2, las=1, tick=TRUE)
        }
    if(i%%nr==0) {
        axis(side=1, at=t_index[(t_index %% 10)==0])
        }
    }
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)

3.4 v2: a unconstrained; b unconstrained

## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL

## number of regr params = level + slope
mm <- 2

## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
  
  ## print loop ID
  print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
  ## get popns
  assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
  ## number of assets
  nn <- dim(assets_mpg)[1]
  
  ## Process eqn
  ## B is Identity
  BB <- "identity"
  ## U is col vec of 0's
  UU <- "zero"
  ## Q will be block diagonal
  QQ <- matrix(list(0),mm*nn,mm*nn)
  ## upper L block for alpha
  aa <- 1:nn
  ## equal var-cov
  # QQ[aa,aa] <- "alpha.cov"
  # diag(QQ[aa,aa]) <- "alpha.var"
  ## diagonal & unequal
  diag(QQ[aa,aa]) <- paste0("alpha",seq(nn))
  for(j in 1:(nn-1)) {
    for(k in (j+1):nn) {
      QQ[j,k] <- paste0("alpha",j,k)
      QQ[k,j] <- paste0("alpha",j,k)
    }
  }
  ## lower R block for beta
  bb <- aa + nn
  ## equal var-cov
  # QQ[bb,bb] <- "beta.cov"
  # diag(QQ[bb,bb]) <- "beta.var"
  ## diagonal & unequal
  diag(QQ[bb,bb]) <- paste0("beta",seq(nn))
  for(j in 1:(nn-1)) {
    for(k in (j+1):nn) {
      QQ[j+nn,k+nn] <- paste0("beta",j,k)
      QQ[k+nn,j+nn] <- paste0("beta",j,k)
    }
  }
  
  ## Observation eqn
  ## Z is array of regr variables (ie, 1's and market returns)
  ZZ <- array(NA, c(nn,mm*nn,TT))
  for(t in 1:TT) {
    ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
  }
  ## A is col vec of 0's
  AA <- "zero"
  ## assume same obs var & cov
  ##RR <- matrix(list(0),nn,nn)
  ##RR[aa,aa] <- "obs.cov"
  ##diag(RR[aa,aa]) <- "obs.var"
  RR <- "diagonal and equal"
  
  ## starting values for regr parameters
  ini_list <- list(x0=matrix(1, mm*nn, 1))
  
  ## list of model matrices & vectors
  mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
  ##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
  
  ## list of control values
  con_list <- list(safe=TRUE, allow.degen=FALSE,
                   conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
  
  ## fit multivariate DLM
  mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
  
  ## get list of Kalman filter output
  kf_out = MARSSkfss(mod_res[[i]])
  
  ## ts of forecasted betas
  betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
  
  ## ts of SE of betas; nxT matrix
  betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
  
  ## ts of alphas
  alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
  ## ts of SE{alphas}
  alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
  
  ## ts of betas
  betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
  ## ts of SE{betas}
  betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
  
}
## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1886 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1886 iterations. 
## Log-likelihood: -350.0825 
## AIC: 810.165   AICc: 837.665   
##  
##           Estimate
## R.diag     0.21763
## Q.alpha1   0.35201
## Q.alpha12  0.28274
## Q.alpha13  0.46744
## Q.alpha14  0.46376
## Q.alpha15  0.33098
## Q.alpha16  0.34731
## Q.alpha2   0.50670
## Q.alpha23  0.40905
## Q.alpha24  0.25186
## Q.alpha25  0.52711
## Q.alpha26  0.32887
## Q.alpha3   0.64873
## Q.alpha34  0.60303
## Q.alpha35  0.49702
## Q.alpha36  0.44837
## Q.alpha4   0.66726
## Q.alpha45  0.31397
## Q.alpha46  0.42392
## Q.alpha5   0.61457
## Q.alpha56  0.38241
## Q.alpha6   0.39532
## Q.beta1    0.08182
## Q.beta12   0.13085
## Q.beta13   0.06558
## Q.beta14   0.03405
## Q.beta15  -0.00181
## Q.beta16   0.04432
## Q.beta2    0.28132
## Q.beta23   0.05297
## Q.beta24   0.03642
## Q.beta25   0.13262
## Q.beta26   0.07967
## Q.beta3    0.09699
## Q.beta34   0.07430
## Q.beta35   0.01698
## Q.beta36   0.07638
## Q.beta4    0.18435
## Q.beta45   0.53115
## Q.beta46   0.24618
## Q.beta5    2.18821
## Q.beta56   0.80106
## Q.beta6    0.34428
## x0.X1      0.56813
## x0.X2      1.97133
## x0.X3      0.82577
## x0.X4      0.25570
## x0.X5      1.44091
## x0.X6      0.64654
## x0.X7      0.91927
## x0.X8      3.62940
## x0.X9     -0.01827
## x0.X10    -0.08166
## x0.X11     1.00034
## x0.X12     0.38951
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 1905 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1905 iterations. 
## Log-likelihood: -344.0416 
## AIC: 798.0832   AICc: 827.0034   
##  
##            Estimate
## R.diag     0.242256
## Q.alpha1   1.408860
## Q.alpha12  1.613053
## Q.alpha13  0.734787
## Q.alpha14  0.674628
## Q.alpha15  1.414349
## Q.alpha16  1.470110
## Q.alpha2   1.874476
## Q.alpha23  0.882608
## Q.alpha24  0.801018
## Q.alpha25  1.615979
## Q.alpha26  1.648587
## Q.alpha3   0.936178
## Q.alpha34  1.026439
## Q.alpha35  0.719197
## Q.alpha36  0.226078
## Q.alpha4   1.165429
## Q.alpha45  0.656512
## Q.alpha46  0.039200
## Q.alpha5   1.420653
## Q.alpha56  1.493419
## Q.alpha6   2.064064
## Q.beta1    0.001881
## Q.beta12  -0.000151
## Q.beta13  -0.001419
## Q.beta14  -0.013128
## Q.beta15  -0.001147
## Q.beta16   0.022776
## Q.beta2    0.002387
## Q.beta23  -0.003008
## Q.beta24   0.005799
## Q.beta25   0.000618
## Q.beta26  -0.006824
## Q.beta3    0.005940
## Q.beta34   0.004748
## Q.beta35   0.000453
## Q.beta36  -0.012415
## Q.beta4    0.103576
## Q.beta45   0.009616
## Q.beta46  -0.173020
## Q.beta5    0.000974
## Q.beta56  -0.015933
## Q.beta6    0.293174
## x0.X1     -0.391740
## x0.X2     -0.870161
## x0.X3     -0.664801
## x0.X4     -0.190172
## x0.X5     -0.380069
## x0.X6     -0.301613
## x0.X7     -0.530841
## x0.X8     -0.243570
## x0.X9     -0.757735
## x0.X10    -0.268795
## x0.X11    -0.692295
## x0.X12    -0.375913
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 1303 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1303 iterations. 
## Log-likelihood: -131.9363 
## AIC: 301.8725   AICc: 308.4242   
##  
##           Estimate
## R.diag       0.121
## Q.alpha1     0.443
## Q.alpha12    0.357
## Q.alpha13    0.386
## Q.alpha2     0.293
## Q.alpha23    0.289
## Q.alpha3     0.428
## Q.beta1      0.192
## Q.beta12     0.153
## Q.beta13     0.142
## Q.beta2      0.122
## Q.beta23     0.113
## Q.beta3      0.105
## x0.X1       -0.582
## x0.X2       -0.581
## x0.X3       -0.246
## x0.X4        0.243
## x0.X5        0.431
## x0.X6        0.357
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1891 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1891 iterations. 
## Log-likelihood: -332.1111 
## AIC: 774.2221   AICc: 801.9699   
##  
##           Estimate
## R.diag     0.24632
## Q.alpha1   1.41757
## Q.alpha12  0.85683
## Q.alpha13  1.02300
## Q.alpha14  1.13028
## Q.alpha15  1.19826
## Q.alpha16  1.75236
## Q.alpha2   0.73748
## Q.alpha23  0.66275
## Q.alpha24  0.71257
## Q.alpha25  0.60989
## Q.alpha26  0.89624
## Q.alpha3   0.74942
## Q.alpha34  0.82637
## Q.alpha35  0.84293
## Q.alpha36  1.25107
## Q.alpha4   0.91557
## Q.alpha45  0.94303
## Q.alpha46  1.41794
## Q.alpha5   1.07329
## Q.alpha56  1.57804
## Q.alpha6   2.46076
## Q.beta1    0.02213
## Q.beta12  -0.01139
## Q.beta13  -0.00460
## Q.beta14   0.01833
## Q.beta15   0.08796
## Q.beta16   0.04825
## Q.beta2    0.00599
## Q.beta23   0.00245
## Q.beta24  -0.00945
## Q.beta25  -0.04558
## Q.beta26  -0.02495
## Q.beta3    0.00104
## Q.beta34  -0.00382
## Q.beta35  -0.01853
## Q.beta36  -0.01011
## Q.beta4    0.01521
## Q.beta45   0.07293
## Q.beta46   0.04002
## Q.beta5    0.35056
## Q.beta56   0.19215
## Q.beta6    0.10543
## x0.X1     -0.67205
## x0.X2     -1.00755
## x0.X3     -0.63267
## x0.X4     -0.55048
## x0.X5     -0.43898
## x0.X6     -0.33243
## x0.X7     -0.57318
## x0.X8     -0.52855
## x0.X9     -0.24796
## x0.X10    -0.79321
## x0.X11    -0.57759
## x0.X12    -1.23466
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop

3.4.1 Plots of results

betas <- betas[,names_mat$ID]

# total num of ts
NN <- dim(betas)[2]

# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
             rainbow(3, start=0.2, end=0.45, v=0.7),
             rainbow(6, start=0.8, end=0.95, v=0.85),
             rainbow(6, start=0.55, end=0.7))

# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7

# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))

# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)

# draw plots
for(i in 1:NN) {
    plot(t_index, betas[,i], type="n",
         ylim=yext, xaxt="n", yaxt="n",
         xlab="", ylab="", main=names_mat[i,"long"])
    rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
    abline(h=0, lty="dashed")
    box()
    lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
    lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
    lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
    #text(xin, yin, names.pop[i], adj=c(0,1))   
    if(i<=nr | nc==1) {
#       axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
#       axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
        axis(side=2, las=1, tick=TRUE)
        }
    if(i%%nr==0) {
        axis(side=1, at=t_index[(t_index %% 10)==0])
        }
    }
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)

3.5 v2: a unconstrained; b equalvarcov

## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL

## number of regr params = level + slope
mm <- 2

## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
  
  ## print loop ID
  print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
  ## get popns
  assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
  ## number of assets
  nn <- dim(assets_mpg)[1]
  
  ## Process eqn
  ## B is Identity
  BB <- "identity"
  ## U is col vec of 0's
  UU <- "zero"
  ## Q will be block diagonal
  QQ <- matrix(list(0),mm*nn,mm*nn)
  ## upper L block for alpha
  aa <- 1:nn
  ## equal var-cov
  # QQ[aa,aa] <- "alpha.cov"
  # diag(QQ[aa,aa]) <- "alpha.var"
  ## diagonal & unequal
  diag(QQ[aa,aa]) <- paste0("alpha",seq(nn))
  for(j in 1:(nn-1)) {
    for(k in (j+1):nn) {
      QQ[j,k] <- paste0("alpha",j,k)
      QQ[k,j] <- paste0("alpha",j,k)
    }
  }
  ## lower R block for beta
  bb <- aa + nn
  ## equal var-cov
  QQ[bb,bb] <- "beta.cov"
  diag(QQ[bb,bb]) <- "beta.var"
  ## diagonal & unequal
  # diag(QQ[bb,bb]) <- paste0("beta",seq(nn))
  # for(j in 1:(nn-1)) {
  #   for(k in (j+1):nn) {
  #     QQ[j+nn,k+nn] <- paste0("beta",j,k)
  #     QQ[k+nn,j+nn] <- paste0("beta",j,k)
  #   }
  # }
  
  ## Observation eqn
  ## Z is array of regr variables (ie, 1's and market returns)
  ZZ <- array(NA, c(nn,mm*nn,TT))
  for(t in 1:TT) {
    ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
  }
  ## A is col vec of 0's
  AA <- "zero"
  ## assume same obs var & cov
  ##RR <- matrix(list(0),nn,nn)
  ##RR[aa,aa] <- "obs.cov"
  ##diag(RR[aa,aa]) <- "obs.var"
  RR <- "diagonal and equal"
  
  ## starting values for regr parameters
  ini_list <- list(x0=matrix(1, mm*nn, 1))
  
  ## list of model matrices & vectors
  mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
  ##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
  
  ## list of control values
  con_list <- list(safe=TRUE, allow.degen=FALSE,
                   conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
  
  ## fit multivariate DLM
  mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
  
  ## get list of Kalman filter output
  kf_out = MARSSkfss(mod_res[[i]])
  
  ## ts of forecasted betas
  betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
  
  ## ts of SE of betas; nxT matrix
  betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
  
  ## ts of alphas
  alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
  ## ts of SE{alphas}
  alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
  
  ## ts of betas
  betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
  ## ts of SE{betas}
  betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
  
}
## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 2003 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2003 iterations. 
## Log-likelihood: -358.099 
## AIC: 788.198   AICc: 799.1609   
##  
##            Estimate
## R.diag       0.2833
## Q.alpha1     0.3294
## Q.alpha12    0.2799
## Q.alpha13    0.4307
## Q.alpha14    0.4216
## Q.alpha15    0.2948
## Q.alpha16    0.3199
## Q.alpha2     0.6465
## Q.alpha23    0.3446
## Q.alpha24    0.1407
## Q.alpha25    0.4108
## Q.alpha26    0.3044
## Q.alpha3     0.6389
## Q.alpha34    0.6086
## Q.alpha35    0.6259
## Q.alpha36    0.5225
## Q.alpha4     0.6851
## Q.alpha45    0.4276
## Q.alpha46    0.4442
## Q.alpha5     1.3891
## Q.alpha56    0.8241
## Q.alpha6     0.5904
## Q.beta.var   0.0675
## Q.beta.cov   0.0675
## x0.X1        0.8072
## x0.X2        0.0792
## x0.X3        1.5557
## x0.X4        1.6881
## x0.X5        1.0960
## x0.X6        0.8313
## x0.X7        0.9286
## x0.X8        1.6950
## x0.X9        1.0510
## x0.X10       1.1087
## x0.X11       0.6499
## x0.X12       1.0357
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 2053 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2053 iterations. 
## Log-likelihood: -348.8623 
## AIC: 769.7246   AICc: 781.2073   
##  
##             Estimate
## R.diag      2.90e-01
## Q.alpha1    1.41e+00
## Q.alpha12   1.56e+00
## Q.alpha13   6.73e-01
## Q.alpha14   6.51e-01
## Q.alpha15   1.43e+00
## Q.alpha16   1.54e+00
## Q.alpha2    1.77e+00
## Q.alpha23   8.30e-01
## Q.alpha24   8.10e-01
## Q.alpha25   1.59e+00
## Q.alpha26   1.63e+00
## Q.alpha3    9.73e-01
## Q.alpha34   1.10e+00
## Q.alpha35   6.56e-01
## Q.alpha36   9.85e-02
## Q.alpha4    1.26e+00
## Q.alpha45   6.27e-01
## Q.alpha46  -5.57e-02
## Q.alpha5    1.46e+00
## Q.alpha56   1.60e+00
## Q.alpha6    2.31e+00
## Q.beta.var  7.07e-05
## Q.beta.cov  3.70e-05
## x0.X1      -5.90e-01
## x0.X2      -9.56e-01
## x0.X3      -2.72e-01
## x0.X4       7.24e-02
## x0.X5      -6.33e-01
## x0.X6      -8.93e-01
## x0.X7      -6.42e-01
## x0.X8      -4.56e-01
## x0.X9      -3.37e-01
## x0.X10     -8.24e-02
## x0.X11     -7.66e-01
## x0.X12     -9.58e-01
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 1120 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1120 iterations. 
## Log-likelihood: -132.3308 
## AIC: 294.6615   AICc: 298.6615   
##  
##            Estimate
## R.diag       0.1204
## Q.alpha1     0.6003
## Q.alpha12    0.4713
## Q.alpha13    0.4842
## Q.alpha2     0.3765
## Q.alpha23    0.3567
## Q.alpha3     0.4764
## Q.beta.var   0.0537
## Q.beta.cov   0.0537
## x0.X1       -0.3677
## x0.X2       -0.4026
## x0.X3        0.0302
## x0.X4        0.3698
## x0.X5        0.6479
## x0.X6        0.6231
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 2106 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2106 iterations. 
## Log-likelihood: -339.5016 
## AIC: 751.0032   AICc: 762.0571   
##  
##             Estimate
## R.diag      3.11e-01
## Q.alpha1    1.45e+00
## Q.alpha12   8.24e-01
## Q.alpha13   1.01e+00
## Q.alpha14   1.19e+00
## Q.alpha15   1.35e+00
## Q.alpha16   1.85e+00
## Q.alpha2    6.55e-01
## Q.alpha23   6.21e-01
## Q.alpha24   6.73e-01
## Q.alpha25   6.47e-01
## Q.alpha26   8.95e-01
## Q.alpha3    7.17e-01
## Q.alpha34   8.32e-01
## Q.alpha35   9.11e-01
## Q.alpha36   1.26e+00
## Q.alpha4    9.85e-01
## Q.alpha45   1.11e+00
## Q.alpha46   1.55e+00
## Q.alpha5    1.33e+00
## Q.alpha56   1.84e+00
## Q.alpha6    2.59e+00
## Q.beta.var  5.28e-05
## Q.beta.cov  2.85e-05
## x0.X1      -5.14e-01
## x0.X2      -1.14e+00
## x0.X3      -6.34e-01
## x0.X4      -3.91e-01
## x0.X5      -7.59e-02
## x0.X6      -4.51e-02
## x0.X7      -4.54e-01
## x0.X8      -6.30e-01
## x0.X9      -2.76e-01
## x0.X10     -6.67e-01
## x0.X11     -2.45e-01
## x0.X12     -9.47e-01
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop

3.5.1 Plots of results

betas <- betas[,names_mat$ID]

# total num of ts
NN <- dim(betas)[2]

# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
             rainbow(3, start=0.2, end=0.45, v=0.7),
             rainbow(6, start=0.8, end=0.95, v=0.85),
             rainbow(6, start=0.55, end=0.7))

# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7

# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))

# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)

# draw plots
for(i in 1:NN) {
    plot(t_index, betas[,i], type="n",
         ylim=yext, xaxt="n", yaxt="n",
         xlab="", ylab="", main=names_mat[i,"long"])
    rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
    abline(h=0, lty="dashed")
    box()
    lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
    lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
    lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
    #text(xin, yin, names.pop[i], adj=c(0,1))   
    if(i<=nr | nc==1) {
#       axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
#       axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
        axis(side=2, las=1, tick=TRUE)
        }
    if(i%%nr==0) {
        axis(side=1, at=t_index[(t_index %% 10)==0])
        }
    }
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)

3.6 v3: a equalvarcov; b unconstrained

## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL

## number of regr params = level + slope
mm <- 2

## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
  
  ## print loop ID
  print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
  ## get popns
  assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
  ## number of assets
  nn <- dim(assets_mpg)[1]
  
  ## Process eqn
  ## B is Identity
  BB <- "identity"
  ## U is col vec of 0's
  UU <- "zero"
  ## Q will be block diagonal
  QQ <- matrix(list(0),mm*nn,mm*nn)
  ## upper L block for alpha
  aa <- 1:nn
  ## equal var-cov
  QQ[aa,aa] <- "alpha.cov"
  diag(QQ[aa,aa]) <- "alpha.var"
  ## lower R block for beta
  bb <- aa + nn
  ## diagonal & unequal
  diag(QQ[bb,bb]) <- paste0("beta",seq(nn))
  for(j in 1:(nn-1)) {
    for(k in (j+1):nn) {
      QQ[j+nn,k+nn] <- paste0("beta",j,k)
      QQ[k+nn,j+nn] <- paste0("beta",j,k)
    }
  }
  
  ## Observation eqn
  ## Z is array of regr variables (ie, 1's and market returns)
  ZZ <- array(NA, c(nn,mm*nn,TT))
  for(t in 1:TT) {
    ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
  }
  ## A is col vec of 0's
  AA <- "zero"
  ## assume same obs var & cov
  ##RR <- matrix(list(0),nn,nn)
  ##RR[aa,aa] <- "obs.cov"
  ##diag(RR[aa,aa]) <- "obs.var"
  RR <- "diagonal and equal"
  
  ## starting values for regr parameters
  ini_list <- list(x0=matrix(1, mm*nn, 1))
  
  ## list of model matrices & vectors
  mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
  ##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
  
  ## list of control values
  con_list <- list(safe=TRUE, allow.degen=FALSE,
                   conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
  
  ## fit multivariate DLM
  mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
  
  ## get list of Kalman filter output
  kf_out = MARSSkfss(mod_res[[i]])
  
  ## ts of forecasted betas
  betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
  
  ## ts of SE of betas; nxT matrix
  betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
  
  ## ts of alphas
  alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
  ## ts of SE{alphas}
  alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
  
  ## ts of betas
  betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
  ## ts of SE{betas}
  betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
  
}
## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1950 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1950 iterations. 
## Log-likelihood: -360.5338 
## AIC: 793.0676   AICc: 804.0305   
##  
##             Estimate
## R.diag        0.4892
## Q.alpha.var   0.3616
## Q.alpha.cov   0.3616
## Q.beta1       0.0382
## Q.beta12      0.0574
## Q.beta13      0.0320
## Q.beta14      0.0248
## Q.beta15      0.0367
## Q.beta16      0.0306
## Q.beta2       0.2277
## Q.beta23     -0.0424
## Q.beta24     -0.1070
## Q.beta25     -0.0641
## Q.beta26     -0.0480
## Q.beta3       0.0845
## Q.beta34      0.1130
## Q.beta35      0.1068
## Q.beta36      0.0857
## Q.beta4       0.1634
## Q.beta45      0.1454
## Q.beta46      0.1158
## Q.beta5       0.1356
## Q.beta56      0.1086
## Q.beta6       0.0871
## x0.X1         0.3421
## x0.X2         0.4242
## x0.X3         0.4651
## x0.X4         0.3915
## x0.X5         0.3153
## x0.X6         0.3377
## x0.X7         0.5544
## x0.X8         1.8386
## x0.X9         0.1294
## x0.X10        0.0441
## x0.X11       -0.0581
## x0.X12        0.2787
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 1923 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1923 iterations. 
## Log-likelihood: -368.2447 
## AIC: 808.4893   AICc: 819.9721   
##  
##             Estimate
## R.diag        0.4578
## Q.alpha.var   1.0297
## Q.alpha.cov   1.0297
## Q.beta1       0.3012
## Q.beta12      0.2121
## Q.beta13     -0.2561
## Q.beta14     -0.1346
## Q.beta15      0.2690
## Q.beta16      0.5721
## Q.beta2       0.1590
## Q.beta23     -0.1796
## Q.beta24     -0.0879
## Q.beta25      0.1911
## Q.beta26      0.3956
## Q.beta3       0.2181
## Q.beta34      0.1151
## Q.beta35     -0.2286
## Q.beta36     -0.4870
## Q.beta4       0.0651
## Q.beta45     -0.1190
## Q.beta46     -0.2607
## Q.beta5       0.2406
## Q.beta56      0.5097
## Q.beta6       1.0917
## x0.X1        -0.1585
## x0.X2        -0.1619
## x0.X3        -0.4460
## x0.X4        -0.2845
## x0.X5        -0.1714
## x0.X6        -0.1760
## x0.X7        -0.2652
## x0.X8         0.2122
## x0.X9        -0.6678
## x0.X10       -0.1822
## x0.X11       -0.3212
## x0.X12       -0.3869
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 1219 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1219 iterations. 
## Log-likelihood: -134.9975 
## AIC: 299.9951   AICc: 303.9951   
##  
##             Estimate
## R.diag        0.1664
## Q.alpha.var   0.2850
## Q.alpha.cov   0.2850
## Q.beta1       0.3401
## Q.beta12      0.2731
## Q.beta13      0.1665
## Q.beta2       0.2209
## Q.beta23      0.1282
## Q.beta3       0.0995
## x0.X1        -0.5229
## x0.X2        -0.5911
## x0.X3        -0.3300
## x0.X4         0.2691
## x0.X5         0.3843
## x0.X6         0.4081
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1934 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1934 iterations. 
## Log-likelihood: -347.5903 
## AIC: 767.1806   AICc: 778.2345   
##  
##             Estimate
## R.diag        0.3472
## Q.alpha.var   0.6749
## Q.alpha.cov   0.6748
## Q.beta1       0.3710
## Q.beta12      0.1475
## Q.beta13      0.2785
## Q.beta14      0.4376
## Q.beta15      0.7304
## Q.beta16      1.1119
## Q.beta2       0.0598
## Q.beta23      0.1144
## Q.beta24      0.1747
## Q.beta25      0.2846
## Q.beta26      0.4466
## Q.beta3       0.2218
## Q.beta34      0.3309
## Q.beta35      0.5285
## Q.beta36      0.8501
## Q.beta4       0.5167
## Q.beta45      0.8579
## Q.beta46      1.3146
## Q.beta5       1.4689
## Q.beta56      2.1651
## Q.beta6       3.3515
## x0.X1        -1.0269
## x0.X2        -1.1166
## x0.X3        -1.0567
## x0.X4        -1.0128
## x0.X5        -1.1300
## x0.X6        -1.1249
## x0.X7        -0.8751
## x0.X8        -0.6575
## x0.X9        -0.4932
## x0.X10       -1.1376
## x0.X11       -1.2226
## x0.X12       -1.9523
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop

3.6.1 Plots of results

betas <- betas[,names_mat$ID]

# total num of ts
NN <- dim(betas)[2]

# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
             rainbow(3, start=0.2, end=0.45, v=0.7),
             rainbow(6, start=0.8, end=0.95, v=0.85),
             rainbow(6, start=0.55, end=0.7))

# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7

# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))

# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)

# draw plots
for(i in 1:NN) {
    plot(t_index, betas[,i], type="n",
         ylim=yext, xaxt="n", yaxt="n",
         xlab="", ylab="", main=names_mat[i,"long"])
    rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
    abline(h=0, lty="dashed")
    box()
    lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
    lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
    lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
    #text(xin, yin, names.pop[i], adj=c(0,1))   
    if(i<=nr | nc==1) {
#       axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
#       axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
        axis(side=2, las=1, tick=TRUE)
        }
    if(i%%nr==0) {
        axis(side=1, at=t_index[(t_index %% 10)==0])
        }
    }
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)

3.7 v3: a equalvarcov; b equalvarcov

## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL

## number of regr params = level + slope
mm <- 2

## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
  
  ## print loop ID
  print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
  ## get popns
  assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
  ## number of assets
  nn <- dim(assets_mpg)[1]
  
  ## Process eqn
  ## B is Identity
  BB <- "identity"
  ## U is col vec of 0's
  UU <- "zero"
  ## Q will be block diagonal
  QQ <- matrix(list(0),mm*nn,mm*nn)
  ## upper L block for alpha
  aa <- 1:nn
  ## equal var-cov
  QQ[aa,aa] <- "alpha.cov"
  diag(QQ[aa,aa]) <- "alpha.var"
  ## lower R block for beta
  bb <- aa + nn
  ## equal var-cov
  QQ[bb,bb] <- "beta.cov"
  diag(QQ[bb,bb]) <- "beta.var"

  ## Observation eqn
  ## Z is array of regr variables (ie, 1's and market returns)
  ZZ <- array(NA, c(nn,mm*nn,TT))
  for(t in 1:TT) {
    ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
  }
  ## A is col vec of 0's
  AA <- "zero"
  ## assume same obs var & cov
  ##RR <- matrix(list(0),nn,nn)
  ##RR[aa,aa] <- "obs.cov"
  ##diag(RR[aa,aa]) <- "obs.var"
  RR <- "diagonal and equal"
  
  ## starting values for regr parameters
  ini_list <- list(x0=matrix(1, mm*nn, 1))
  
  ## list of model matrices & vectors
  mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
  ##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
  
  ## list of control values
  con_list <- list(safe=TRUE, allow.degen=FALSE,
                   conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
  
  ## fit multivariate DLM
  mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
  
  ## get list of Kalman filter output
  kf_out = MARSSkfss(mod_res[[i]])
  
  ## ts of forecasted betas
  betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
  
  ## ts of SE of betas; nxT matrix
  betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
  
  ## ts of alphas
  alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
  ## ts of SE{alphas}
  alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
  
  ## ts of betas
  betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
  ## ts of SE{betas}
  betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
  
}
## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1916 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1916 iterations. 
## Log-likelihood: -366.8833 
## AIC: 767.7665   AICc: 770.1024   
##  
##             Estimate
## R.diag        0.5887
## Q.alpha.var   0.3221
## Q.alpha.cov   0.3221
## Q.beta.var    0.0502
## Q.beta.cov    0.0501
## x0.X1         0.4183
## x0.X2         0.4331
## x0.X3         0.5491
## x0.X4         0.4894
## x0.X5         0.3806
## x0.X6         0.4096
## x0.X7         0.4309
## x0.X8         0.6753
## x0.X9         0.5201
## x0.X10        0.8094
## x0.X11        0.3756
## x0.X12        0.7383
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 1937 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1937 iterations. 
## Log-likelihood: -376.9972 
## AIC: 787.9943   AICc: 790.4326   
##  
##              Estimate
## R.diag       6.26e-01
## Q.alpha.var  1.07e+00
## Q.alpha.cov  1.07e+00
## Q.beta.var   1.10e-04
## Q.beta.cov   6.04e-05
## x0.X1       -5.73e-01
## x0.X2       -5.90e-01
## x0.X3       -7.16e-01
## x0.X4       -6.62e-01
## x0.X5       -5.82e-01
## x0.X6       -6.32e-01
## x0.X7       -5.60e-01
## x0.X8       -4.34e-01
## x0.X9       -7.52e-01
## x0.X10      -6.55e-01
## x0.X11      -6.77e-01
## x0.X12      -5.62e-01
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 1618 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1618 iterations. 
## Log-likelihood: -138.2772 
## AIC: 298.5543   AICc: 300.6833   
##  
##             Estimate
## R.diag         0.210
## Q.alpha.var    0.338
## Q.alpha.cov    0.338
## Q.beta.var     0.097
## Q.beta.cov     0.097
## x0.X1         -0.494
## x0.X2         -0.534
## x0.X3         -0.407
## x0.X4          0.189
## x0.X5          0.444
## x0.X6          0.433
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1992 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1992 iterations. 
## Log-likelihood: -359.999 
## AIC: 753.998   AICc: 756.3519   
##  
##              Estimate
## R.diag       4.97e-01
## Q.alpha.var  1.01e+00
## Q.alpha.cov  1.01e+00
## Q.beta.var   8.00e-05
## Q.beta.cov   4.74e-05
## x0.X1       -4.03e-01
## x0.X2       -5.21e-01
## x0.X3       -5.00e-01
## x0.X4       -3.95e-01
## x0.X5       -4.40e-01
## x0.X6       -5.11e-01
## x0.X7       -4.63e-01
## x0.X8       -5.31e-01
## x0.X9       -2.76e-01
## x0.X10      -6.60e-01
## x0.X11      -3.44e-01
## x0.X12      -7.99e-01
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop

3.7.1 Plots of results

betas <- betas[,names_mat$ID]

# total num of ts
NN <- dim(betas)[2]

# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
             rainbow(3, start=0.2, end=0.45, v=0.7),
             rainbow(6, start=0.8, end=0.95, v=0.85),
             rainbow(6, start=0.55, end=0.7))

# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7

# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))

# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)

# draw plots
for(i in 1:NN) {
    plot(t_index, betas[,i], type="n",
         ylim=yext, xaxt="n", yaxt="n",
         xlab="", ylab="", main=names_mat[i,"long"])
    rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
    abline(h=0, lty="dashed")
    box()
    lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
    lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
    lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
    #text(xin, yin, names.pop[i], adj=c(0,1))   
    if(i<=nr | nc==1) {
#       axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
#       axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
        axis(side=2, las=1, tick=TRUE)
        }
    if(i%%nr==0) {
        axis(side=1, at=t_index[(t_index %% 10)==0])
        }
    }
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)

3.8 v4: a diag & unequal; b diag & unequal

## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL

## number of regr params = level + slope
mm <- 2

## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
  
  ## print loop ID
  print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
  ## get popns
  assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
  ## number of assets
  nn <- dim(assets_mpg)[1]
  
  ## Process eqn
  ## B is Identity
  BB <- "identity"
  ## U is col vec of 0's
  UU <- "zero"
  ## Q will be block diagonal
  QQ <- matrix(list(0),mm*nn,mm*nn)
  ## upper L block for alpha
  aa <- 1:nn
  ## equal var-cov
  # QQ[aa,aa] <- "alpha.cov"
  # diag(QQ[aa,aa]) <- "alpha.var"
  ## diagonal & unequal
  diag(QQ[aa,aa]) <- paste0("alpha",seq(nn))
  ## lower R block for beta
  bb <- aa + nn
  ## diagonal & unequal
  diag(QQ[bb,bb]) <- paste0("beta",seq(nn))

  ## Observation eqn
  ## Z is array of regr variables (ie, 1's and market returns)
  ZZ <- array(NA, c(nn,mm*nn,TT))
  for(t in 1:TT) {
    ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
  }
  ## A is col vec of 0's
  AA <- "zero"
  ## assume same obs var & cov
  ##RR <- matrix(list(0),nn,nn)
  ##RR[aa,aa] <- "obs.cov"
  ##diag(RR[aa,aa]) <- "obs.var"
  RR <- "diagonal and equal"
  
  ## Initial states
  ##x0 <- rbind(matrix("pi.a",nn,1), matrix("pi.b",nn,1))
  ##V0 <- 100*diag(mm*nn)
  ##x0 <- rbind(matrix(0,nn,1), matrix(1,nn,1))
  ##V0 <- matrix(list(0),mm*nn,mm*nn)
  ##diag(V0[aa,aa]) <- "var.a"
  ##diag(V0[bb,bb]) <- "var.b"
  
  ## starting values for regr parameters
  ini_list <- list(x0=matrix(1, mm*nn, 1))
  
  ## list of model matrices & vectors
  mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
  ##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
  
  ## list of control values
  con_list <- list(safe=TRUE, allow.degen=FALSE,
                   conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
  
  ## fit multivariate DLM
  mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
  
  ## get list of Kalman filter output
  kf_out = MARSSkfss(mod_res[[i]])
  
  ## ts of forecasted betas
  betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
  
  ## ts of SE of betas; nxT matrix
  betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
  
  ## ts of alphas
  alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
  ## ts of SE{alphas}
  alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
  
  ## ts of betas
  betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
  ## ts of SE{betas}
  betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
  
}
## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1378 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1378 iterations. 
## Log-likelihood: -418.6651 
## AIC: 887.3302   AICc: 892.4483   
##  
##           Estimate
## R.diag    8.43e-01
## Q.alpha1  6.52e-05
## Q.alpha2  1.32e-01
## Q.alpha3  2.11e-02
## Q.alpha4  1.01e-01
## Q.alpha5  8.29e-05
## Q.alpha6  5.10e-02
## Q.beta1   2.41e-02
## Q.beta2   1.59e-01
## Q.beta3   1.26e-01
## Q.beta4   2.72e-04
## Q.beta5   4.27e-02
## Q.beta6   7.89e-03
## x0.X1    -6.42e-02
## x0.X2     1.07e+00
## x0.X3    -1.95e-01
## x0.X4    -1.10e-01
## x0.X5    -8.11e-02
## x0.X6    -1.90e-01
## x0.X7    -4.61e-02
## x0.X8     2.25e+00
## x0.X9    -1.07e+00
## x0.X10   -1.11e-01
## x0.X11   -1.62e-01
## x0.X12   -4.89e-02
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 3158 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 3158 iterations. 
## Log-likelihood: -458.1922 
## AIC: 966.3844   AICc: 971.7341   
##  
##           Estimate
## R.diag    1.45e+00
## Q.alpha1  1.16e-02
## Q.alpha2  4.30e-01
## Q.alpha3  1.04e-04
## Q.alpha4  1.23e-01
## Q.alpha5  4.24e-05
## Q.alpha6  4.23e-05
## Q.beta1   1.07e-04
## Q.beta2   1.01e-04
## Q.beta3   6.94e-05
## Q.beta4   1.11e-04
## Q.beta5   2.75e-02
## Q.beta6   7.91e-02
## x0.X1    -4.03e-01
## x0.X2    -1.31e+00
## x0.X3    -2.40e-01
## x0.X4    -6.99e-01
## x0.X5    -1.67e-01
## x0.X6    -3.34e-01
## x0.X7    -3.17e-01
## x0.X8    -3.57e-01
## x0.X9    -4.62e-01
## x0.X10   -5.31e-01
## x0.X11   -1.90e-01
## x0.X12   -2.84e-01
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 1259 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1259 iterations. 
## Log-likelihood: -172.4242 
## AIC: 370.8483   AICc: 373.832   
##  
##           Estimate
## R.diag    4.27e-01
## Q.alpha1  1.08e-02
## Q.alpha2  1.59e-03
## Q.alpha3  7.53e-05
## Q.beta1   3.87e-01
## Q.beta2   2.22e-01
## Q.beta3   1.25e-01
## x0.X1    -3.68e-01
## x0.X2    -1.41e-01
## x0.X3     1.64e-01
## x0.X4     6.99e-02
## x0.X5     9.89e-01
## x0.X6     7.69e-01
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## 
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1331 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1331 iterations. 
## Log-likelihood: -462.4741 
## AIC: 974.9483   AICc: 980.107   
##  
##           Estimate
## R.diag    3.92e-01
## Q.alpha1  1.06e+00
## Q.alpha2  5.17e-01
## Q.alpha3  5.60e-01
## Q.alpha4  5.79e-01
## Q.alpha5  5.83e-01
## Q.alpha6  2.42e+00
## Q.beta1   1.53e-04
## Q.beta2   9.44e-05
## Q.beta3   9.76e-05
## Q.beta4   9.17e-05
## Q.beta5   5.42e-01
## Q.beta6   3.28e-04
## x0.X1    -5.04e-01
## x0.X2    -1.05e+00
## x0.X3    -8.87e-01
## x0.X4    -2.71e-01
## x0.X5    -4.24e-02
## x0.X6     6.98e-02
## x0.X7    -4.12e-01
## x0.X8    -5.86e-01
## x0.X9    -3.01e-01
## x0.X10   -4.91e-01
## x0.X11   -1.33e-01
## x0.X12   -8.43e-01
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.
## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop

3.8.1 Plots of results

betas <- betas[,names_mat$ID]

# total num of ts
NN <- dim(betas)[2]

# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
             rainbow(3, start=0.2, end=0.45, v=0.7),
             rainbow(6, start=0.8, end=0.95, v=0.85),
             rainbow(6, start=0.55, end=0.7))

# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7

# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))

# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
          ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)

# draw plots
for(i in 1:NN) {
    plot(t_index, betas[,i], type="n",
         ylim=yext, xaxt="n", yaxt="n",
         xlab="", ylab="", main=names_mat[i,"long"])
    rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
    abline(h=0, lty="dashed")
    box()
    lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
    lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
    lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
    #text(xin, yin, names.pop[i], adj=c(0,1))   
    if(i<=nr | nc==1) {
#       axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
#       axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
        axis(side=2, las=1, tick=TRUE)
        }
    if(i%%nr==0) {
        axis(side=1, at=t_index[(t_index %% 10)==0])
        }
    }
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)